MASSACHUSETTS INSTITUTE OF TECHNOLOGY 
ARTIFICIAL INTELLIGENCE LABORATORY 



A.I. Memo 989. November, 1987 



Rigidity and smoothness of motion. 



S. Ullman and A. YuiHe 



Abstract. Many theories of structure from motion divide the process into 
two parts which are solved using different assumptions. Smoothness of the 
velocity field is often assumed to solve the motion correspondence problem 
and then rigidity is typically used to recover the 3D structure. We prove 
results showing that, in a statistical sense, smoothness of the velocity field 
foUows from rigidity of the motion. 

(c) Massachusetts Institute of Technology, 1987 

This report describes research done at the Artificial Intelligence Labo- 
ratory of the Massachusetts Institute of Technology. Support for the lab- 
oratory's artificial intelhgence research is provided in part by the Advanced 
Research Projects Agency of the Department of Defense under Office of Naval 
Research contract N00014-85-K-0124. 



Yuille Ullman 

RIGIDITY AND SMOOTHNESS OF MOTION 



The smoothness assumption in measuring visual motion 

The problem of measuring visual motion is, in many situations, undercon- 
strained. That is, the information in the changing image is insufficent to 
determine the motion uniquely. This indeterminateness is often referred to as 
the "aperture problem" (for example Man and Ullman, 1982). An additional 
constraint is therefore required for resolving the ambiguity and determining 
the velocity field uniquely. An important constraint that has been proposed 
for solving the problem is a smoothness constraint. (Fennema and Thomp- 
son 1979, Horn and Schunck 1981, Hildreth 1984, Nakayama and Silverman 
1986). When applied to the problem of measuring the velocity of image con- 
tours, this smoothness assumption has been formulated by Hildreth in the 
foUowing manner. Of aU the velocity fields that are consistent with the trans- 
forming image, select the velocity field that minimizes the overall variation. 
That is, if v denotes the velocity at a point, and s the arc length along the 
curve, th« preferred velocity field is the one that minimizes the integral: 

along the curve (the principle of least variation). It has been shown that 
this method of determining the velocity field gives good results under a wide 
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range of conditions, and it seems to correspond in many Ccises to the velocity 
field perceived by human observers (Hildreth 1986, Nakayama and Silverman 
1986). 

Justifying the smoothness assumption. 

The main rationale raised in support of the smoothness assumption is that 

the velocity field induced by smooth contours in motion is expected to be 

smooth. This argument is insufficent, however, for justifying the principle of 

least variation. The smoothness assumption implies that if the object is known 

to be smooth, the velocity field should not contain discontinuities. It does not 

imply, however, that "the smoother the better". Consider, for instance, two 

points in the image separated by a small distance, and moving with image 

velocities Vi and V2. One may argue that if they lie on the same surface, 

extremely high values of H^i - F2II are unlikely, because of some physical 

limitations on the motion of objects in space. But why should one assume 

that a relative velocity of, say, 0.2 deg/sec is less likely than 0.1 deg/sec? Such 

a preference is assumed by the least variation principle (that favors velocity 

fields in which the difference in velocity between neighboring points is as small 

as possible), but cannot be defended only on the basis of a general smoothness 

argument. Clearly, a stronger constraint than just the general smoothness of 
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surfaces is required. 

Another general property that seems to provide a useful constraint in the 
analysis of visual motion is rigidity. Computational studies have shown that 
the 3D stucture of rigid and quasi-rigid objects can be recovered by looking 
for the most rigid interpretation possible of the changing image (Ullman 1979, 
1984). 

In this paper we argue that local rigidity of the object and the principle of 
least variation in the velocity field are related. To investigate this we consider 
a rod moving in three-dimensional space. The rod is allowed to move in a semi- 
rigid manner: it can rotate and translate freely and is also allowed a certain 
expansion. We calculate the relative velocity between the two endpoints of 
the rod, as projected on the image plane. The calculations show that under 
a wide range of conditions this velocity distribution is peaked at zero velocity 
and decreases monotonically at higher velocities. As a result, the projected 
velocity field of points linked rigidly together is Ukely to be consistent with 
the principle of least variation in the velocity field. There are two factors 
contributing to this strong bias towards small differential velocities in the 
projected velocity field: (i) the rigidity of the link between the two points and 
(ii) the effect of projection from 3-D to the 2-D image plane. 

If we imagine an object being made up of a collection of semi-rigid rods 
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this analysis suggests that the projected velocity field is likely to be smooth. 
We conclude that the principle of least variation is a reasonable constraint 
that can be justified for the projected velocity field of rigid or locally rigid 
objects. 

Section 1. The two-dimensional case. 

First we consider a rigid rod moving in two dimensions and being pro- 
jected onto a line. Because it is rigid the rod's velocity can be split into 
rotational and translational components. For our purposes the translational 
component can be ignored as it will not contribute to the derivative of the 
velocity field on the line. 

Let the projected length of the rod be R and the real length be r, the 
angle to the vertical be and the angular velocity be u. Then the projected 
velocity distribution $(u) is 

$(u) = J 6{R- rsine)6(u - njcosd)^r{r)^U^)drdujde (1.1) 
where S denotes the Dirac delta function, $^(r) and $^(a;) are the distribution 
of r and w respectively. rsinO and rucos9 are the projected lengths and 
velocity respectively. We assume that rods are equally likely to have any 
length between and r^^x and set $r(r) = k between and Vmax (with 
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We can integrate with respect to using the first delta function. This 
gives (see Appendix 1) 

, ^ _ ^^ ^iu - uy/r^ - R^)^Uu)drdu;. (1.2) 
We integrate with respect to r to obtain 

$(w) = / z^^J. _ n Q\ 

where w^^^ = u/^/rl~^^R?. This lower bound comes from the delta func- 
tion integration. 

Now we examine the behavior of $(w). If we differentiate it with respect 
to u we get 

So, /or any probability distribution $u;(w) (which by definition must be a 
positive function) the two terms on the right hand side of (1.4) must be 
negative. $(«) must decrease strictly monotonically with u and has a unique 
maximum at w = 0. 

Thus whatever the distribution of the rotation the most likely projected 
velocity is zero. This result is similar to that obtained by Ullman (1979) 
for the motion of dots. It was shown that if the motion of a random dot is 
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described by an isotropic probability distribution function in three dimensions 
then the probability distribution for the projected two-dimensional motion is 
peaked at 0. 

Note that our result is independent of whether we choose an upper bound 
for r. This can easily be seen by setting Vmax ^ oo and u; k^ in (1.4). In 
fact having an upper bound for r makes d^/du more negative and makes 
$(w) decay faster. Henceforth we assume for simplicity that there is no upper 
bound for r. 

To see the connection with the smoothness assumption observe that in 
the limit Hildreth's smoothness measure can be expressed as 

/ll£lP''-E9^ (1.5) 

where the sum can be taken over a set of rigid rotating rods. The results above 
show that for each rod considered independently the probability distribution 
of the u is peaked at zero. Thus, with this independence assumption, we 
argue that the most likely distribution of the whole contour is the one that 
minimizes (1.5). We wiU relax the independence assumption in section 4. 

Section 2. The three-dimensional case. 

We now extend the analysis to a rod moving in 3-space. Consider a rod 
of length r projected into the image plane, which has unit normal vector k. 
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The rod's direction in 3-D space is denoted by r which is a unit vector in the 
direction r. The rod is projected into a vector fp where 

Tp = f- {f-k)k. (2.1) 

The rod's rotation is described by a vector a;, where u is the frequency and 
tD is the axis of rotation. The velocity v is given by 

v = cDxr. (2.2) 

The projected velocity field Vp is given by 

Vp = v-{v-k)k. (2.3) 

We consider rods with projected length R and projected velocity tj. We are 
interested in the projected velocity distribution $(f/). Since the projected 
velocity distribution is rotationally symmetric we can express this as U^{U) 
(see Appendix 1). If the rods' length distribution is $^(f) and their rotation 
distribution is $u;(^) then the projected velocity distribution is given by 

$(w) ^ j KR- \rp\)KU - \vp\)^r{i^^u;(Cj)dfdu}. (2.4) 

Care is needed in specifying the domain of integration of (2.4). From (2.2) 
we see that only the component of u; perpendicular to f contributes to the 
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velocity. We must require w • f = 0. We impose this by substituting a delta 
function into (2.4) and then integrating over all uj. 
Thus we have 

U^U) = JhR- Kmij 'fjSiU- \vj,\)(^^(Cj)dfduj. (2.5) 

We assume all directions of the rod are equally likely and the rods are 
also equally likely to have any length. So we set 

^(r) = 1. (2.6) 

We also assume that the rotation is equally likely to be in any direction. Then 
$u;(^) is a function of |u;| = lj only. 

Um) = Js{R- \fp\)6{u . J^6(U - \vp\)^^{u;)dfdu. (2.7) 

Choose an orthogonal basis I J, k. Define angles 0, (p, a, /3. These angles 
are illustrated in figure 1. ^ is the angle between the unit vector and the k 
axis, (fi is the angle between the unit vector projected onto the x,y plane and 
the X-axis. Similarly for a and /3. It follows that 



f = sin$ cos(pi + sind sirup] + cosdk 

(2.8) 
a; = sina cos(3i + sina sin (3 J + cosak 
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(i) 




(ii) 



k i 




Figure 1 Definitiona of the angles, see text. 



Then 



du = du sina dad/S 
df= dr 3in0 dOdip 



and 



(2.9) 



|fj,| = rainO. 
Define the angle between w and f to be r by 



(2.10) 



COST =z uj ' r. 



(2.11) 
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Then from (2.8) we have 
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COST = cos9 cosa -f sinff sinot cos((p - f3). (2.12) 

Define the angle p by 






From (2.11) we obtain 



cospsinr = (tD x f • ib). (2.14) 

From (2.8) we find 

cospsinr = sin0sina{(p — /?) (2.15) 

and we calculate 

\vp\ = ujsiuT s'm p. (2.16) 

We use (2.10), (2.11) and (2.16) to write the delta functions as 

H^-\rp\)S{oj'r)6{U-\vj,\) = 6{R-r sine)6(cosT)6{U -QrsiriT sinp). (2.17) 
We now need to change the variables of integration. We have 
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dujdf=: dijjdr sinasinS dadj3d0d(p. (2.18) 

The integrand only depends on (f and jS as (p- f3. We can set ^ = y? - /? and 
integrate out /3. Then we have 

COST = CO30 cosa 4- sinB sinacostp 

cospsinr = sin$ sina sinip • (2.19) 

dQdr = -wdudr sinasinO dadij^dO 
Now we must change the variables of integration from a, V to r, p: 



da da \ 
dr dp \ 

dadip = det^ ^ I drdp. (2.20) 



^ Z 






From (2.19) we eliminate the terms containing ^ to obtain 

cosoL = COST cosO ± sina(sin^$ - cos^p)K (2.21) 

Define functions F(T,p,a,ip) andC?(r,^,a, V') by 



(2.22) 



^{'''yP^ot^t/)) = COST - cosBcosa - sinO sinacosip 

^(^» pyOi,il;) = cosp siuT - sinO sina sini/) 
It is a standard result of Jacobian transformations (it can be obtained by 

using formula (2.20) twice, once changing variables from a, ^ to F,G and 

then changing F, G to r, p) that 
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dadtj) = det 



From (2.22) we obtain 



dF dF 



dG dG 



/det 



dF dF 
dG dG 



drdp. 



(2.23) 



and 



det 



dF 
■57 



dG dG 
"W Tp 



— Sin T Sin p 



(2.24) 



det 



dF dF 



dG 



dG 



Thus we calculate 



= sinO sina{sin0 cosa - cosfp cosB sinot) . (2.25) 



dadxl) = 



sm^T stnp 



sinB sina{sin0 cosa — cost/) cosd sina} ^ ^* 
We substitute for cos^ from (2.19) and obtain 



(2.26) 



dad'tj; = 



then using (2.21) we find 



sm T sinp 



sina{cosa — cosOcosr} 



drdp 



dadfp = 



sinr stnp 



sina(sin^9 — cos^ p)h 



drdp. 



Hence 



(2.27) 



(2.28) 



• /I J j; sinr sinp sin9 
stna stnddadtl) = drdp. 

(sin^B — cos'^p)h 



(2.29) 
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We substitute (2.17), (2.18) and (2.29) into (2.7) 



U^{u) = / 6(R - rsin0)S{cosT)6(u - ursinrsinp) 

... sinr sinpsinO , , ,^ 

^'^^^^ 7~~T^ — xi ^^ '^Z' ^^ ^ «^^ 
(^in-^^ — cos^pP 

Now we do the integration with respect to r to obtain 



(2.30) 



TTih(TT\ f ^(^ ~ '''sinB)6(U — (jjrsinp) , ^ ^ . 

Um) = J ^,i„,g_\„,2p)i '^' 'iU^)stnpsmedpd»d^dr. (2.31) 

We integrate with respect to 9 to obtain 

U^(u)= f ^^^ ~ '^''^'''P)^<^i^)^iMR/r)dp(L;dr 

J {r^-R'')h{B?/T^-cos^p)h ' ^^-^^^ 

We integrate with respect to p and obtain 



^{IJ)= { R^u;{i^)(LjdT 



(2.33) 



r(r2 - i22)i(^2^^2 + ^2 _ a;2r2)i(a,2r2 _ ^2)^ ' 
The domain of integration is restricted to the region of the a; - r plane where 
the integrand is real (for simplicity we have omitted the bounds of integration 
in the previous equations). It is specified by 

r>R 
R^iJ + «^ > J^t"^ (2.34) 

a;r > u. 
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Change the variable of integration from r to X where 
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iJr'^ = w2 ^. X 



(2.35) 



then 



with domain of integration 



X < R^Lj^ 



X>0 



(2.37) 



w' + X > u;2i22 



Define a new variable Y by 



y2 = u;^r2 _ X 



(2.38) 



then 



^u) = y 



Ru^^{u;)du;dY 



{iJ^R^ - y2)^(^2 _ y2)^(^2 ^. ^2 J22 _ y2) 



-r^ (2.39) 



with 
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y >0 

Y <u;R (2.40) 

Y <u. 
We can further simplify this expression by the substitutions Y = usinO 

and u = {u/R)Q. This yields 






^e^^'ia^-sin^ey/^in^ + cos^ey (2.41) 

It is clear from this expression that $(w) will decrease with u if and only 
if ^u;M does not grow faster than uj^. Thus for all realistic cases $(w) will 
be a monotonic decreasing function of u. 

Section 3. Expansion and Contraction. 

We now consider the situation in two dimensions when the rod is allowed 
to expand and contract as well as rotate. We can write 



r = rr 



(3.1) 



and differentiating gives 



r = rr-\- rOh (32) 

where n is the unit vector perpendicular to f . We define the expansion coef- 
ficient s by 
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r = rs 
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(3.3) 



and write (3.2) as 

f = rsr + r0h. (3.4^ 

The projected velocity is then given by 

U = rssinS + rujcosB (3.5) 

where is the angle that f makes with the normal to the image plane and 
= u. The projected velocity distribution is given by 

HU) = J S{R - rsin0)6{U - (r s sinO + ru;cos$))^,(3)^^{uj) dr duj ds dS 

(3.6) 
where, as in Section 1, we have assumed that $(r) = 1 so the rods have no 
preferred lengths. We do the integral of ^U) with respect to to obtain 

m)=J(r'-R')-'^'SiU-sR-u;(r'-Ry/^)^^(s)^^(^)drdujds. (3.7) 
We now integrate with respect to r to find 

m) = j{{U - sRf + uj^R'')-'f''^M^^{^)d^ds. (3.8) 



YuiUe Ullman 

Care must be taken to ensure that this integral is evaluated over the correct 
limits, those for which the integrand is well defined. For a; > we need s < ^ 
and for u; < we have s > ^. This domain of integration is shown in figure 
(2). We can split the integral up into two parts 



HU) = ^i{U)-\-^2{U) 



Figure 2 The domain of iniegraiion, see text. 
where - - 



(3.9) 




ru/soo f»^U/R 



^i(U) = / du {(U- sRf + B?iJr^'''^.(s)^^(u)ds (3. 



10) 



and 
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MU) = / duj {{U- sRf + ie2u;2)-i/2$,(5)$,(a;)rf^. (3.II) 



u;=-oo Js=U/R 

We impose the conditions 



$5(5) = ^,{-s) (3.12) 

and 

*u;H = M-'*')- (3.13) 

These ensure that contraction is as likely as expansion and that rotation is 

equally likely to be any direction. 

We use (3.12) and (3.13) to rewrite (3.9) as 



^U) = du {(U- sRf + R}uj'')-^'H,{3)^^{uj)ds 

J Lj=Q ./»= — 00 

f(jJ=oo rs=-U/R 

+ / ^^ {iP + sRf + i22u;2)-i/2$,(5)$^(a;)rf5. (3.14^) 

./u/=0 Ja=— 00 

We now substitute t = Rs/u and a; = uQ/R and obtain 



•fl=oo rt=l 



Note that if we set ^s(ut/R) = 6(ut/R) we recover (1.3). 
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We differentiate (3.14a) to obtain 



^^ J u=0 Js=—oo 

ru)=oo fS= — U/R 

- duj {{U - sRf +B?u'^)-^f\U -sR)^,{s)^^(uj)ds. (3.15) 

./a;=0 Ja=—(x> 



We can simplify this expression by setting s to 5 - 2C//i2 in the second term 
on the right hand side. This yields 

^$ y.a<=oo ta=-U/R ^ 

dU=- Lo "^""Js^-oo iiU + sR)2 + R2^2y/2 i^R+U)i^s{s)-^s{s+2U/R))^^(^)ds. 

(3.16) 

In the range of integration (u + sR) is always negative. So the integral will 
be always negative provided $,(5) - $,(5 + 2u/R) is negative for s < -u/R. 
If we make 5 positive, using (3.12), this becomes 

^s(s - 2U/R) - $,(5) > 0, (3.17) 

which is true for all monotonically decreasing functions. Then we have |^ < 
everywhere for any function $a;(a;). Hence the projected velocity distribution 
will once more be peaked about U = 0. 

K we relax the conditions (3.17) we are no longer assured of having a 
monotonically decreasing projected velocity distribution for all $^(a;). If 
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$,(5) is linear in s we obtain a flat distribution 

mi) = const. (3.18) 

If $5(5) is quadratic in 5 we find 



du 



= -2jd^^^ + JJdu;dss\s'+u;')-y^^^^i^) (3.19) 
which will only be negative for some distributions $w ((<;). 

Thus allowing the rod to expand and contract reduces the strength of the 
result of Section 1. But if we make the reasonable assumption that $,(5) is 
at most linear in s the result will still hold. 

Section 4. Curves from rods. 

The situation becomes more complex when we consider several rods 
joined together to form a curve. The results above show that for each rod 
individually the projected velocity distribution is peaked at zero. If the rods 
are joined, however, these probability distributions are no longer independent. 
We now sketch an argument suggesting that this loss of independence does 
not affect the result. 

Suppose we have a set of probability distributions Pi(xi), i = 1,..N 
which are all peaked at zero. We now define the joint probability distribution 
P{x) where x = (ari,X2, ...,xn). If we impose no constraints this distribution 
is given by 
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P{x)dx = IliPi(xi)Uidxi. (4.1) 

Now impose a consistency constraint J^j ^t = 0. We now have 

J P{x)dx = j YiiPi{xi)6{Y, Xi)nidxi. (4.2) 

The x.'s can no longer vary independently and the delta function imposes the 
consistency constraint. We can get integrate out the delta function thereby 
reducing the number of independent variable io N - I. It is convenient to 
choose these to be Xi - x^+i, for i = 1,...,7V - 1. We can now write (4.2) as 

j P{x)dx = J UiPiixi - Xi+i )Ilid{xi - Xi+i ). (4.3) 

We can set y,- = Xi - Xj+i, i = 1, ...,iV - 1 and rewrite 

Piif) = Myi)P2(y2)"-PN-i{yN-i)PN{-yi - 2/2 - ... - vn-i). (4.4) 

It is straightforward to see that the maximum value of P(y) occurs when the 
Vi are all zero. Thus P is still peaked at zero even with the constraints. 

Section 5. An alternative approach. 

The results derived in the last few sections support the assumption (Hil- 
dreth 1984) that the velocity field along a contour is smooth. This assumpti 
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is used to solve the aperture problem by prescribing the smoothest possible 
velocity field consistent with the data. 

An alternative approach to the aperture problem was taken by Waxman 
and Wohn (1986). They showed that if the object was locally planar and 
moving rigidly then, assuming perspective projection, the local velocity field 
on the image plane would be locally quadratic in the x,y coordinates. They 
then defined regions in which these expansions were valid and checked for 
consistency between these regions. In these regions they could do a least 
squares best fit to find the coefficients of these quadratic polonomials, and 
hence the velocity field. 

This approach, however, can also be used to justify the smoothest velocity 
field assumption. The orthographic projection of the motion of a rigid body 
will obey 



Vx- A + By-^ Cz, (5.1a) 



Vy = D-\-Ex^- Fz, (5.16) 

where A, B,C,D,E,F are constants and z is the (unknown) depth coordinate. 
If we also assume the object is planar then z obeys 
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z = px -^ qy + s, (5.2) 

where p,q,r are constants. Substituting (5.2) into (5.1) gives linear expres- 
sions for Vx and Vy in terms of x and y. So if we assume local planarity and 
local rigidity (in the style of Waxman and Wohn (1986)) the velocity fields 
will locally be linear in x and y. The measure of smoothest velocity, /(v), 
used by Hildreth is the integral along the curve of the function 

where 5 is the arc length. The velocity fields will be locally linear if and only 
if J{v) is zero. Thus we can think of the smoothest velocity field approach as 
a local method of assuming local rigidity and local planarity. 

Section 6. Conclusion. 

The arguments in the first four sections of this paper suggest that locally 
rigid, or semi-rigid objects, will tend to project a smooth velocity field i 
the image. Moreover, assuming random motions and limited expansion 
contraction, this field will tend to be as smooth as possible. In the fifth 
section we noted that rigidity of an object and local smoothness of its surface 
wiU also lead to a smooth image motion. 



m 



or 
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These arguments support the view that maximixing smoothness is a good 
heuristic to use for motion correspondence and that it is a sensible way to solve 
the aperture problem. 
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Appendix 

We first describe the general method for integrating delta functions. 
Suppose we have an integral /(a, 6) 

I{a,b)= I f(x)6{x - xo)dx (A.l) 

J a 

where 6{x) is the Dirac delta function, f{x) is an arbitrary function, xq an 
arbitrary point and h>a. The value of the integral is 

I{0'M = f{xQ\ if xoe[a,b] {A.2a) 

I{a,b) = 0, otherwise. (^4.26) 

This result can be generalized to integrals of form 
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J(a,6) = / f{x)S(g(x) - c)dx {A.3) 

J a 

where g{x) is an arbitrary function and c an arbitrary number. All points Xi 
with g{xi) = c will contribute to this interval. Suppose there is only one such 
point. If there are several we can divide the integral up into regions with only 
one such point. Consider one such point ar = 0. The function g{x) can be 
expanded in a Taylor series about this point 

g(x) = c + g'{Qi){x) + 0{x'^). {AA) 

If we change the coordinate to u where u — g'(0)x we can write the 
integral as 

J(a, 6) = / f{u/g'iO))S{u + 0(u'))-^. (A.5) 

Jg'iO)a g'{0) 

The value of the integral will depend on the sign of g'(Q). If it is negative 
the bounds of the integral will be reversed and the integral will change sign. 
Therefore 



J(a,6) = /(0)j~^, 0€[a,6] {A.Qa) 



J(a, b) = 0, otherwise. (AM) 
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We now consider the form of the probability distribution function for a 
rotational symmetric vector in two-dimensional space. Suppose the function 
is ^Q{Cj)duj. If it is rotationally symmetric (and hence depends only on the 
modulus a; of u5) it can be written 

^ci{Cj)doj = $a;(w)wc?a>rfy>, (^-7) 

by changing to radial coordinates a;, (fi in the lj space. The (p component can 
be integrated out to give a distribution LJ^u{u)d(jj. 
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